library(raster)
library(data.table)
library(dplyr)
library(plm)


path<-"XXXX/HYDE Files/"

# desired<-c("1100AD","1200AD","1300AD","1400AD","1500AD","1600AD","1700AD","1800AD","1900AD","2000AD")
# NB - The Netherlands Environmental Assessment Agency has changed the url for the HYDE data since the first draft of the paper. Below is the current url
# for (d in 1:length(desired)){
#   temp<-tempfile()
#   source_html<-"https://dataportaal.pbl.nl/downloads/HYDE/HYDE3.2/2016_beta_release.zip"
#   # source_html<-paste(source_html,desired[d],sep="")
#   # source_html<-paste(source_html,"pop.zip",sep="_")
#   download.file(source_html,temp,method="libcurl")
# 
#   wanted_file<-paste("popd",desired[d],sep="_")
#   wanted_file<-paste(wanted_file,"asc",sep=".")
# 
#   temp_file<-unzip(temp, wanted_file)
#   unlink(temp)
# 
# }

temp<-paste(path,list.files(path)[1],sep="")

temp<-raster(temp)

temp<-aggregate(temp,fac=12,fun=mean,na.rm=T)

projection(temp)<-"+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"

coords<-round(coordinates(temp)-.5)

tsdata<-matrix(NA,dim(coords)[1],length(list.files(path)))

for (i in 1:length(list.files(path))){
  file<-paste(path,list.files(path)[i],sep="")
  popd<-raster(file)
  popd<-aggregate(popd,fac=12,fun=mean,na.rm=T)
  tsdata[,i]<-values(popd)
  
        
}

tsdata<-data.frame(coords,tsdata)

names(tsdata)<-c("Longitude","Latitude","Pop1100AD","Pop1200AD","Pop1300AD","Pop1400AD","Pop1500AD","Pop1600AD","Pop1700AD","Pop1800AD","Pop1900AD","Pop2000AD")


###Import HCED and locating it within cells and centuries

hced<-read.csv("XXXX/HCED Data.csv")
hced<-hced[!is.na(hced$Latitude),]

x<-rasterToPolygons(temp)

swx<-round(coordinates(x)[,1]-.5)
swy<-round(coordinates(x)[,2]-.5)
ID<-paste(swx,swy,sep="-")

x@data[,2]<-ID
identifier<-as.numeric(as.factor(ID))
x@data[,3]<-identifier
names(x@data)[2]<-"ID"
names(x@data)[3]<-"Identifier"


###Converting the coordinates of the battles to a spatial points dataframe

points<-SpatialPointsDataFrame(cbind(hced$Longitude,hced$Latitude),data=hced,proj4string=CRS("+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"))

###Placing the battles in polygons

out<-over(points,x)

hced_ID<-out[,2]

###Placing the battles in centuries

centuries<-seq(1100,2000,100)
start<-centuries[1:length(centuries)-1]
end<-centuries[2:length(centuries)]
hced$Year<-as.numeric(hced$Year)
century<-matrix(NA,length(hced$Year),1)
century<-ifelse(dplyr::between(hced$Year,start[1],end[1])==TRUE,"12th Century",NA)

for (i in 1:length(start)){
  num<-i+1
  cen<-paste("1",as.character(num),sep="")
  cen<-paste(cen,"th Century")
  print(cen)
  century<-ifelse(dplyr::between(hced$Year,start[i],end[i])==TRUE,cen,century)
  
}

century<-ifelse(century=="110 th Century","20th Century",century)

cells<-unique(hced_ID)

###Creating a number of battles variable for each century-cell

cens<-na.omit(unique(century))

for (i in 1:length(cens)){
a<-apply(as.matrix(cells),1,function(x) length(na.omit(hced$ID[hced_ID==x&century==cens[i]])))
assign(paste("battles",cens[i],sep="-"),a)
}

tsdata$ID<-paste(tsdata$Longitude,tsdata$Latitude,sep="-")

bats<-data.frame(cells,`battles-12 th Century`,`battles-13 th Century`,`battles-14 th Century`,`battles-15 th Century`,`battles-16 th Century`,`battles-17 th Century`,`battles-18 th Century`,`battles-19 th Century`,`battles-20th Century`)
names(bats)<-c("ID","bat12","bat13","bat14","bat15","bat16",'bat17',"bat18","bat19","bat20")

tsdata<-left_join(tsdata,bats,by="ID")

nordhaus<-read.csv("XXXXX/Gecon40_post_final .csv")
nord<-data.frame(nordhaus$LAT,nordhaus$LONGITUDE,nordhaus$COUNTRY)
rm(nordhaus)
names(nord)<-c("LAT","LONGITUDE","COUNTRY")
nord$ID<-paste(nord$LONGITUDE,nord$LAT,sep="-")

tsdata<-left_join(tsdata,nord,by="ID")
tsdata<-distinct(tsdata,ID,.keep_all=TRUE)

x<-merge(x,tsdata,by.x="ID",by.y="ID")

region<-matrix(NA,dim(x@data)[1],1)
region[which(x@data$COUNTRY=="American Samoa")]<-"Asia Pacific"
region[which(x@data$COUNTRY=='South Korea')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Philippines')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Australia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Laos')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Samoa')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Brunei Darussalam')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Macao SAR, China')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Singapore')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Cambodia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Malaysia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Solomon Islands')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='China')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Marshall Islands')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Taiwan')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Fiji')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Micronesia, Fed. Sts.')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Marshall Islands')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Thailand')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='French Polynesia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Mongolia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Timor Leste')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Guam')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Myanmar')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Papua New Guinea')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Hong Kong SAR, China')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Nauru')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Tonga')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Indonesia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='New Caledonia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Tuvalu')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Japan')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='New Zealand')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Vanuatu')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Kiribati')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Northern Mariana Is.')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Vietnam')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='North Korea')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Palau')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Albania')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Gibraltar')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Norway')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Andorra')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Greece')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Poland')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Armenia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Greenland')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Portugal')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Austria')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Hungary')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Romania')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Azerbaijan')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Iceland')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Russia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Belarus')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Ireland')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='San Marino')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Belgium')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Isle of Man')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Serbia and Montenegro')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Bosnia&Herzegovina')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Italy')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Slovakia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Bulgaria')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Kazakhstan')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Slovenia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Kosovo')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Spain')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Croatia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Kyrgyztan')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Sweden')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Cyprus')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Latvia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Switzerland')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Czech Republic')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Liechtenstein')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Tajikistan')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Denmark')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Lithuania')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Turkey')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Estonia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Luxembourg')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Turkmenistan')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Moldova')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Ukraine')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Finland')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Monaco')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='United Kingdom')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='France')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Montenegro')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Uzbekistan')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Georgia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Netherlands')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Germany')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='North Macedonia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Antigua and Barbuda')]<-'South America'
region[which(x@data$COUNTRY=='Curacao')]<-'South America'
region[which(x@data$COUNTRY=='Paraguay')]<-'South America'
region[which(x@data$COUNTRY=='Argentina')]<-'South America'
region[which(x@data$COUNTRY=='Dominica')]<-'South America'
region[which(x@data$COUNTRY=='Peru')]<-'South America'
region[which(x@data$COUNTRY=='Aruba')]<-'South America'
region[which(x@data$COUNTRY=='Dominican Republic')]<-'South America'
region[which(x@data$COUNTRY=='Puerto Rico')]<-'South America'
region[which(x@data$COUNTRY=='Bahamas')]<-'South America'
region[which(x@data$COUNTRY=='Ecuador')]<-'South America'
region[which(x@data$COUNTRY=='Sint Maarten')]<-'South America'
region[which(x@data$COUNTRY=='Barbados')]<-'South America'
region[which(x@data$COUNTRY=='El Salvador')]<-'South America'
region[which(x@data$COUNTRY=='St. Kitts and Nevis')]<-'South America'
region[which(x@data$COUNTRY=='Belize')]<-'South America'
region[which(x@data$COUNTRY=='Grenada')]<-'South America'
region[which(x@data$COUNTRY=='Guadeloupe')]<-'South America'
region[which(x@data$COUNTRY=='St. Lucia')]<-'South America'
region[which(x@data$COUNTRY=='Bolivia')]<-'South America'
region[which(x@data$COUNTRY=='Guatemala')]<-'South America'
region[which(x@data$COUNTRY=='St. Martin')]<-'South America'
region[which(x@data$COUNTRY=='Brazil')]<-'South America'
region[which(x@data$COUNTRY=='Guyana')]<-'South America'
region[which(x@data$COUNTRY=='St. Vincent and the Grenadines')]<-'South America'
region[which(x@data$COUNTRY=='British Virgin Islands')]<-'South America'
region[which(x@data$COUNTRY=='Haiti')]<-'South America'
region[which(x@data$COUNTRY=='Suriname')]<-'South America'
region[which(x@data$COUNTRY=='Cayman Islands')]<-'South America'
region[which(x@data$COUNTRY=='Honduras')]<-'South America'
region[which(x@data$COUNTRY=='Trinidad and Tobago')]<-'South America'
region[which(x@data$COUNTRY=='Chile')]<-'South America'
region[which(x@data$COUNTRY=='Jamaica')]<-'South America'
region[which(x@data$COUNTRY=='Turks and Caicos Islands')]<-'South America'
region[which(x@data$COUNTRY=='Colombia')]<-'South America'
region[which(x@data$COUNTRY=='Mexico')]<-'South America'
region[which(x@data$COUNTRY=='Uruguay')]<-'South America'
region[which(x@data$COUNTRY=='Costa Rica')]<-'South America'
region[which(x@data$COUNTRY=='Nicaragua')]<-'South America'
region[which(x@data$COUNTRY=='Venezuela')]<-'South America'
region[which(x@data$COUNTRY=='Cuba')]<-'South America'
region[which(x@data$COUNTRY=='Panama')]<-'South America'
region[which(x@data$COUNTRY=='Virgin Islands (U.S.)')]<-'South America'
region[which(x@data$COUNTRY=='Algeria')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Jordan')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Qatar')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Bahrain')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Kuwait')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Saudi Arabia')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Djibouti')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Lebanon')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Syria')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Egypt')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Libya')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Tunisia')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Iran')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Malta')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='United Arab Emirates')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Iraq')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Morocco')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='West Bank and Gaza')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Israel')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Oman')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Yemen')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Bermuda')]<-'North America'
region[which(x@data$COUNTRY=='Canada')]<-'North America'
region[which(x@data$COUNTRY=='United States')]<-'North America'
region[which(x@data$COUNTRY=='Afghanistan')]<-'South Asia'
region[which(x@data$COUNTRY=='India')]<-'South Asia'
region[which(x@data$COUNTRY=='Pakistan')]<-'South Asia'
region[which(x@data$COUNTRY=='Bangladesh')]<-'South Asia'
region[which(x@data$COUNTRY=='Maldives')]<-'South Asia'
region[which(x@data$COUNTRY=='Sri Lanka')]<-'South Asia'
region[which(x@data$COUNTRY=='Ile de Reunion')]<-'South Asia'
region[which(x@data$COUNTRY=='Bhutan')]<-'South Asia'
region[which(x@data$COUNTRY=='Nepal')]<-'South Asia'
region[which(x@data$COUNTRY=='Angola')]<-'Africa'
region[which(x@data$COUNTRY=='Ethiopia')]<-'Africa'
region[which(x@data$COUNTRY=='Niger')]<-'Africa'
region[which(x@data$COUNTRY=='Benin')]<-'Africa'
region[which(x@data$COUNTRY=='Gabon')]<-'Africa'
region[which(x@data$COUNTRY=='Nigeria')]<-'Africa'
region[which(x@data$COUNTRY=='Botswana')]<-'Africa'
region[which(x@data$COUNTRY=='Gambia')]<-'Africa'
region[which(x@data$COUNTRY=='Rwanda')]<-'Africa'
region[which(x@data$COUNTRY=='Burkina Faso')]<-'Africa'
region[which(x@data$COUNTRY=='Ghana')]<-'Africa'
region[which(x@data$COUNTRY=='Sao Tome and Principe')]<-'Africa'
region[which(x@data$COUNTRY=='Burundi')]<-'Africa'
region[which(x@data$COUNTRY=='Guinea')]<-'Africa'
region[which(x@data$COUNTRY=='Senegal')]<-'Africa'
region[which(x@data$COUNTRY=='Cabo Verde')]<-'Africa'
region[which(x@data$COUNTRY=='Guinea Bissau')]<-'Africa'
region[which(x@data$COUNTRY=='Seychelles')]<-'Africa'
region[which(x@data$COUNTRY=='Cameroon')]<-'Africa'
region[which(x@data$COUNTRY=='Kenya')]<-'Africa'
region[which(x@data$COUNTRY=='Sierra Leone')]<-'Africa'
region[which(x@data$COUNTRY=='Central African Republic')]<-'Africa'
region[which(x@data$COUNTRY=='Lesotho')]<-'Africa'
region[which(x@data$COUNTRY=='Somalia')]<-'Africa'
region[which(x@data$COUNTRY=='Chad')]<-'Africa'
region[which(x@data$COUNTRY=='Liberia')]<-'Africa'
region[which(x@data$COUNTRY=='South Africa')]<-'Africa'
region[which(x@data$COUNTRY=='Comoros')]<-'Africa'
region[which(x@data$COUNTRY=='Madagascar')]<-'Africa'
region[which(x@data$COUNTRY=='South Sudan')]<-'Africa'
region[which(x@data$COUNTRY=='Democratic Republic of Congo')]<-'Africa'
region[which(x@data$COUNTRY=='Malawi')]<-'Africa'
region[which(x@data$COUNTRY=='Sudan')]<-'Africa'
region[which(x@data$COUNTRY=='Congo')]<-'Africa'
region[which(x@data$COUNTRY=='Mali')]<-'Africa'
region[which(x@data$COUNTRY=='Tanzania')]<-'Africa'
region[which(x@data$COUNTRY=="Cote d'Ivoire")]<-'Africa'
region[which(x@data$COUNTRY=='Mauritania')]<-'Africa'
region[which(x@data$COUNTRY=='Togo')]<-'Africa'
region[which(x@data$COUNTRY=='Equatorial Guinea')]<-'Africa'
region[which(x@data$COUNTRY=='Mauritius')]<-'Africa'
region[which(x@data$COUNTRY=='Uganda')]<-'Africa'
region[which(x@data$COUNTRY=='Eritrea')]<-'Africa'
region[which(x@data$COUNTRY=='Mozambique')]<-'Africa'
region[which(x@data$COUNTRY=='Zambia')]<-'Africa'
region[which(x@data$COUNTRY=='Eswatini')]<-'Africa'
region[which(x@data$COUNTRY=='Namibia')]<-'Africa'
region[which(x@data$COUNTRY=='Zimbabwe')]<-'Africa'



x@data[,28]<-region
names(x@data)[28]<-"Region"

x@data$Region[which(x@data$COUNTRY=='Uzbekistan')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Tajikistan')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Kazakhstan')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Turkmenistan')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Kyrgyztan')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Armenia')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Georgia')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Azerbaijan')]<-'Central Asia'


library(spdep)

###Creating the neighbors' list

x.temp <- poly2nb(x)
nb2INLA("temp.graph", x.temp)
W.adj <- paste(getwd(),"/temp.graph",sep="")

###Creating the time series dataset

ID<-rep(x@data$ID,times=1,each=9)
Identifier<-rep(x@data$Identifier,times=1,each=9)
region<-rep(x@data$Region,times=1,each=9)
year<-rep(seq(1200,2000,100),times=length(x@data$Identifier),each=1)

data<-data.frame(ID,Identifier,year,region)

battles<-matrix(NA,dim(data)[1],1)

battles[data$year==1200]<-x@data$bat12
battles[data$year==1300]<-x@data$bat13
battles[data$year==1400]<-x@data$bat14
battles[data$year==1500]<-x@data$bat15
battles[data$year==1600]<-x@data$bat16
battles[data$year==1700]<-x@data$bat17
battles[data$year==1800]<-x@data$bat18
battles[data$year==1900]<-x@data$bat19
battles[data$year==1200]<-x@data$bat20

popd<-matrix(NA,dim(data)[1],1)

popd[data$year==1200]<-x@data$Pop1200AD
popd[data$year==1300]<-x@data$Pop1300AD
popd[data$year==1400]<-x@data$Pop1400AD
popd[data$year==1500]<-x@data$Pop1500AD
popd[data$year==1600]<-x@data$Pop1600AD
popd[data$year==1700]<-x@data$Pop1700AD
popd[data$year==1800]<-x@data$Pop1800AD
popd[data$year==1900]<-x@data$Pop1900AD
popd[data$year==2000]<-x@data$Pop2000AD

data<-data.frame(data,battles,popd)

data<-pdata.frame(data,index=c("Identifier","year"))

data$battles[is.na(data$battles)]<-0

data$battles_diff<-data$battles-lag(data$battles,100)

data$pop_diff<-data$popd-lag(data$popd,100)

data$pop_diff1<-lag(data$pop_diff,100)
data$pop_diff2<-lag(data$pop_diff,200)
data$pop_diff3<-lag(data$pop_diff,300)

data$battles_diff1<-lag(data$battles_diff,100)
data$battles_diff2<-lag(data$battles_diff,200)
data$battles_diff3<-lag(data$battles_diff,300)

data$Identifier<-as.numeric(data$Identifier)

write.csv(data,file="Time Series Data.csv")

library(INLA)

###Model 1###

f.inla <- battles_diff ~ pop_diff1 + f(Identifier, model="bym", graph=W.adj) 
set.seed(1234)
mod <- inla(f.inla,family="gaussian",data=data)


mean<-summary(mod)[[3]]["pop_diff1",1]
sd<-summary(mod)[[3]]["pop_diff1",2]
sim_beta<-rnorm(1000,mean,sd)
png("Global Model 1.png")
plot(density(sim_beta),xlab=expression(beta),main="Effect of Lagged Population Density Growth on Battles")
segments(summary(mod)[[3]]["pop_diff1",3],0,summary(mod)[[3]]["pop_diff1",3],max(density(sim_beta)$y),lty=2)
segments(summary(mod)[[3]]["pop_diff1",5],0,summary(mod)[[3]]["pop_diff1",5],max(density(sim_beta)$y),lty=2)
dev.off()


###Model 2###


f.inla <- pop_diff ~ battles_diff1 + f(Identifier, model="bym", graph=W.adj) 
set.seed(1234)
mod <- inla(f.inla,family="gaussian",data=data)


mean<-summary(mod)[[3]]["battles_diff1",1]
sd<-summary(mod)[[3]]["battles_diff1",2]
sim_beta<-rnorm(1000,mean,sd)
png("Global Model 2.png")
plot(density(sim_beta),xlab=expression(beta),main="Effect of Lagged Battles on Population Density Growth")
segments(summary(mod)[[3]]["battles_diff1",3],0,summary(mod)[[3]]["battles_diff1",3],max(density(sim_beta)$y),lty=2)
segments(summary(mod)[[3]]["battles_diff1",5],0,summary(mod)[[3]]["battles_diff1",5],max(density(sim_beta)$y),lty=2)
dev.off()

###Model 3###


f.inla <- pop_diff ~ battles_diff1 + pop_diff2 + f(Identifier, model="bym", graph=W.adj) 
set.seed(1234)
mod <- inla(f.inla,family="gaussian",data=data)


mean<-summary(mod)[[3]]["battles_diff1",1]
sd<-summary(mod)[[3]]["battles_diff1",2]
sim_beta<-rnorm(1000,mean,sd)
png("Global Model 3.png")
plot(density(sim_beta),xlab=expression(beta),main="Effect of Lagged Battles on Population Density Growth")
segments(summary(mod)[[3]]["battles_diff1",3],0,summary(mod)[[3]]["battles_diff1",3],max(density(sim_beta)$y),lty=2)
segments(summary(mod)[[3]]["battles_diff1",5],0,summary(mod)[[3]]["battles_diff1",5],max(density(sim_beta)$y),lty=2)
dev.off()

###By Region

regions<-levels(as.factor(data$region))
means<-NULL;lowers<-NULL; uppers<-NULL

for (i in 1:length(regions)){

reg<-regions[i]
dat<-subset(data,data$region==reg)
loc <- unique(dat$ID)
temp_ID<- 1:length(loc)
new_dat<-data.frame(loc,temp_ID)
dat <- merge(dat,new_dat,by.x="ID",by.y="loc")
dat <- dat[order(dat$temp_ID),]


x.temp <- x[x$ID%in%sort(unique(dat$ID)),]
x.new<-x.temp[order(x.temp$Identifier),]
x.new <- poly2nb(x.new)
nb2INLA("temp_new.graph", x.new)
W.adj_new <- paste(getwd(),"/temp_new.graph",sep="")
f.inla <- pop_diff ~ battles_diff2 + pop_diff3 + f(temp_ID, model="bym", graph=W.adj_new) 
set.seed(1234)
mod <- inla(f.inla,family="gaussian",data=dat)


mean<-summary(mod)[[3]]["battles_diff2",1]
sd<-summary(mod)[[3]]["battles_diff2",2]
sim_beta<-rnorm(1000,mean,sd)

mean<-summary(mod)[[3]][2,1]
lower<-summary(mod)[[3]][2,3]
upper<-summary(mod)[[3]][2,5]

means<-c(means,mean)
lowers<-c(lowers,lower)
uppers<-c(uppers,upper)

title<-gsub("/"," ",reg)
title<-paste(title,"Model .png",sep=" ")
png(title)
plot(density(sim_beta),xlab=expression(beta),main="Effect of Lagged Battles on Population Density Growth")
segments(summary(mod)[[3]]["battles_diff2",3],0,summary(mod)[[3]]["battles_diff2",3],max(density(sim_beta)$y),lty=2)
segments(summary(mod)[[3]]["battles_diff2",5],0,summary(mod)[[3]]["battles_diff2",5],max(density(sim_beta)$y),lty=2)
dev.off()

}


out<-data.frame(regions,means,lowers,uppers)






